Comme pour toute entreprise de location, il est important de connaître la demande que le service peut rencontrer. Dans un but évident d'optimisation des coûts. Dans notre exemple, une société de location de vélos, il peut être intéressant de prédire le nombre de vélos loués par heure. Notamment pour anticiper les stocks, le personnel en service ou encore positionner les vélos dans des endroits de la ville connaissant une forte demande. C'est dans cette optique, que nous allons chercher à comprendre les features qui impactent le nombre de vélos loués par heure.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
file = pd.read_csv("data.csv")
file
print(file.info(),"\n")
print(file.isna().sum())
import seaborn as sns
import plotly.express as px
import plotly.graph_objects as go
fig, ax = plt.subplots(figsize=(20,15))
matrix = np.triu(file.corr())
cmap = sns.diverging_palette(220, 20, sep=20, as_cmap=True)
sns.heatmap(file.corr(),annot = True, fmt='.1g', mask=matrix, ax=ax, cmap=cmap)
We can see few strong correlations, between the cnt and if a customer is registered or casual. So, they can overfit the model.
In addition, features temp and atemp have a correlation of 1. So, we can keep only one feature to represent both in our model. This will made a model more faster to compute.
Seasons and month have obviously a hight correlation.
fig = px.violin(file, x="weekday", y="casual", box=True)
fig.show()
from scipy import stats
print('Average :',file.groupby(file.weekday)['casual'].mean())
print('Standard Deviation :',file.groupby(file.weekday)['casual'].std())
fig = px.violin(file, x="weekday", y="registered", box=True)
fig.show()
print('Average :',file.groupby(file.weekday)['registered'].mean())
print('Standard Deviation :',file.groupby(file.weekday)['registered'].std())
Si nous prenons le jour 0 (lundi) du cas "registered". Avec une moyenne de 121 locations/lundi et une std de 106, cela signifie que nous avons 68% de chance que le nbr de location soit entre (121-106)-(121+106) = 15-207.
Les autres jours suivent une tendance similaire. Ce qui est un écart assez important. On peut donc penser que le dataset possède un certain nombre de "outliers" et que donc la médiane est une donnée plus précise que la moyenne dans notre cas.
#Dropping the outlier rows with standard deviation. Trimmed mean can be a solution too
factor = 3
upper_lim = file['temp'].mean () + file['temp'].std () * factor
lower_lim = file['temp'].mean () - file['temp'].std () * factor
file = file[(file['temp'] < upper_lim) & (file['temp'] > lower_lim)]
upper_lim = file['atemp'].mean () + file['atemp'].std () * factor
lower_lim = file['atemp'].mean () - file['atemp'].std () * factor
file = file[(file['atemp'] < upper_lim) & (file['atemp'] > lower_lim)]
upper_lim = file['hum'].mean () + file['hum'].std () * factor
lower_lim = file['hum'].mean () - file['hum'].std () * factor
file = file[(file['hum'] < upper_lim) & (file['hum'] > lower_lim)]
upper_lim = file['windspeed'].mean () + file['windspeed'].std () * factor
lower_lim = file['windspeed'].mean () - file['windspeed'].std () * factor
file = file[(file['windspeed'] < upper_lim) & (file['windspeed'] > lower_lim)]
upper_lim = file['registered'].mean () + file['registered'].std () * factor
lower_lim = file['registered'].mean () - file['registered'].std () * factor
file = file[(file['registered'] < upper_lim) & (file['registered'] > lower_lim)]
upper_lim = file['casual'].mean () + file['casual'].std () * factor
lower_lim = file['casual'].mean () - file['casual'].std () * factor
file = file[(file['casual'] < upper_lim) & (file['casual'] > lower_lim)]
file = file.reset_index()
file = file.drop(['index'], axis=1)
file
day = file.groupby(['weekday']).sum()
day = pd.DataFrame(day)
day = day.reset_index()
fig = go.Figure()
fig.update_layout(
title="Distribution de la location des vélos par jour",
xaxis_title="Jour",
yaxis_title="Nbr de vélos",
)
fig.add_trace(go.Scatter(
x=day['weekday'],
y=day['cnt'],
))
fig.show()
On peut voir une augmentation graduelle de l'utilisation du service au cours de la semaine. Avec un pic le samedi qui pourrait correspondre à une utilisation plus importantes des utilisateurs "casual"
trip = file.groupby(file.hr)['cnt'].sum()
trip = pd.DataFrame(trip)
trip = trip.reset_index()
fig = go.Figure()
fig.update_layout(
title="Distribution de la location des vélos par heure",
xaxis_title="Heure",
yaxis_title="Nbr de vélos",
)
fig.add_trace(go.Scatter(
x=trip['hr'],
y=trip['cnt'],
))
fig.show()
En premier lieu, on peut observer une utilisation du service évolutive en fonction de l'heure de la journée. Avec notamment, 2 fortes hausses d'utilisation à 8H et 19H. Ce qui correspond aux horaires de début et fin de travail.
trip = file
trip = trip[trip['workingday'] == 0]
trip = trip.groupby(trip.hr)['cnt'].sum()
trip2 = file
trip2 = trip2[trip2['workingday'] == 1]
trip2 = trip2.groupby(trip2.hr)['cnt'].sum()
trip = pd.DataFrame(trip)
trip2 = pd.DataFrame(trip2)
trip = trip.reset_index()
trip2 = trip2.reset_index()
fig = go.Figure()
fig.update_layout(
title="Distribution de la location des vélos par heure",
xaxis_title="Heure",
yaxis_title="Nbr de vélos",
)
fig.add_trace(go.Scatter(
x=trip['hr'],
y=trip['cnt'],
name="Week-End"
))
fig.add_trace(go.Scatter(
x=trip2['hr'],
y=trip2['cnt'],
name="Working Day"
))
fig.show()
Le WE, on peut voir une utilisation du service assez uniforme entre 10H et 20H. Alors que les jours de travail la courbe suit la même tendance que le plot juste au dessus. Avec des tendances moins smooth et plus affirmées. Dans le premier plot, la baisse entre 9H et 16H est moins importante car compensée par les chiffres du WE, qui connaissent une augmentation "uniforme" là où la semaine connait une baisse.
Dans la prédiction de vélo, il sera donc important de se situer si nous sommes dans un "working day" ou non.
casu = file.groupby(file.weekday)['casual'].sum()
regi = file.groupby(file.weekday)['registered'].sum()
casu = pd.DataFrame(casu)
casu = casu.reset_index()
regi = pd.DataFrame(regi)
regi = regi.reset_index()
fig = go.Figure()
fig.add_trace(go.Bar(
x=regi['weekday'],
y=regi['registered'],
name='Registered Clients',
marker_color='#12c2e9'
))
fig.add_trace(go.Bar(
x=casu['weekday'],
y=casu['casual'],
name='Casual Clients',
marker_color='#c471ed'
))
fig.update_layout(barmode='group', xaxis_tickangle=-45, title_text="Distribution of features 'Casual' & 'Registered' per day")
fig.show()
casu = file.groupby(file.mnth)['casual'].sum()
regi = file.groupby(file.mnth)['registered'].sum()
casu = pd.DataFrame(casu)
casu = casu.reset_index()
regi = pd.DataFrame(regi)
regi = regi.reset_index()
fig = go.Figure()
fig.add_trace(go.Bar(
x=regi['mnth'],
y=regi['registered'],
name='Registered Clients',
marker_color='#12c2e9'
))
fig.add_trace(go.Bar(
x=casu['mnth'],
y=casu['casual'],
name='Casual Clients',
marker_color='#c471ed'
))
fig.update_layout(barmode='group', xaxis_tickangle=-45, title_text="Distribution of features 'Casual' & 'Registered' per month")
fig.show()
casu = file.groupby(file.season)['casual'].sum()
regi = file.groupby(file.season)['registered'].sum()
casu = pd.DataFrame(casu)
casu = casu.reset_index()
regi = pd.DataFrame(regi)
regi = regi.reset_index()
fig = go.Figure()
fig.add_trace(go.Bar(
x=regi['season'],
y=regi['registered'],
name='Registered Clients',
marker_color='#12c2e9'
))
fig.add_trace(go.Bar(
x=casu['season'],
y=casu['casual'],
name='Casual Clients',
marker_color='#c471ed'
))
fig.update_layout(barmode='group', xaxis_tickangle=-45, title_text="Distribution of features 'Casual' & 'Registered' per season")
fig.show()
casu = file.groupby(file.yr)['casual'].sum()
regi = file.groupby(file.yr)['registered'].sum()
casu = pd.DataFrame(casu)
casu = casu.reset_index()
regi = pd.DataFrame(regi)
regi = regi.reset_index()
fig = go.Figure()
fig.add_trace(go.Bar(
x=regi['yr'],
y=regi['registered'],
name='Registered Clients',
marker_color='#12c2e9'
))
fig.add_trace(go.Bar(
x=casu['yr'],
y=casu['casual'],
name='Casual Clients',
marker_color='#c471ed'
))
fig.update_layout(barmode='group', xaxis_tickangle=-45, title_text="Distribution of features 'Casual' & 'Registered' per year")
fig.show()
Les utilisateurs sont très majoritairements des personnes abonnées au service. Le Lundi et Dimanche ainsi que les saisons/mois chauds connaissent une hausse des clients non abonnés, mais cela reste moindre.
La distribution par jour, sauf pour le lundi est dimanche, est assez proche.
Le service semble avoir du succès car d'une année à une autre, le nombre de location a augmenté de 1.2M à 2M. Soit une augmentation de 66.67%
day = file.groupby(['hum']).sum()
day = pd.DataFrame(day)
day = day.reset_index()
fig = go.Figure()
fig.update_layout(
title="Distribution de la location des vélos en fonction de l'humidité",
xaxis_title="Humidité",
yaxis_title="Nbr de vélos",
)
fig.add_trace(go.Scatter(
x=day['hum'],
y=day['cnt'],
))
fig.show()
On ne peut observer de pattern concernant l'influence de l'humidité avec la nombre de location
day = file.groupby(['windspeed']).sum()
day = pd.DataFrame(day)
day = day.reset_index()
fig = go.Figure()
fig.update_layout(
title="Distribution de la location des vélos en fonction de la force du vent",
xaxis_title="Windspeed",
yaxis_title="Nbr de vélos",
)
fig.add_trace(go.Scatter(
x=day['windspeed'],
y=day['cnt'],
))
fig.show()
On peut voir que l'ulisation du service diminue plus le vent devient important. On peut expliquer la baisse entre 0 et 0.9 par le fait qu'il est rare d'avoir aucun ou très peu de vent.
fig = px.sunburst(file, path=['weekday','weathersit'], values='cnt')
fig.show()
Peut importe le jour de la semaine, les utilisateurs n'utilisent pas le service au-delà du niveau 2 de météo. C'est à dire quand il pleut.
from sklearn.cluster import KMeans
from sklearn import cluster
from sklearn.preprocessing import StandardScaler
temp = file.groupby(file.temp)['cnt'].sum()
temp = pd.DataFrame(temp)
temp = temp.reset_index()
fig = px.bar(temp, x='temp', y='cnt',color='cnt',height=400)
fig.show()
x = temp['temp']
y = temp['cnt']
temp = StandardScaler().fit_transform(temp)
kmeans = KMeans(n_clusters=3)
kmeans.fit(temp)
cluster = kmeans.predict(temp)
plt.scatter(x,y, c=cluster)
On peut voir que le cluster vert, qui correspond aux données entre 0.30 et 0.8 sur le diagramme. Représente une hausse de l'utilisation du service pour les températures comprises entre :
mnth = pd.DataFrame(file)
mask = (mnth['dteday'] > '2012-01-01') & (mnth['dteday']<= '2012-12-31')
mnth = mnth.loc[mask]
mnth = mnth.groupby(mnth.mnth)['cnt'].sum()
mnth = mnth.reset_index()
fig = px.bar(mnth, x='mnth', y='cnt',color='cnt',height=400)
fig.show()
from scipy.stats import normaltest
from scipy.stats import shapiro
stat, p = normaltest(mnth['cnt'])
print('Statistics=%.2f, p=%.2f' % (stat, p))
# interpret
alpha = 0.05
if p > alpha:
print('Sample looks Gaussian (fail to reject H0)')
else:
print('Sample does not look Gaussian (reject H0)')
stat, p = shapiro(mnth['cnt'])
print('Statistics=%.3f, p=%.3f' % (stat, p))
# interpret
alpha = 0.05
if p > alpha:
print('Sample looks Gaussian (fail to reject H0)')
else:
print('Sample does not look Gaussian (reject H0)')
#Z-score
mean = mnth['cnt'].mean()
std = mnth['cnt'].std()
print("mean = ",mean)
print("std = ",std)
z = (x - mean)/std
z1 = (107656-mean)/std
z2 = (147656-mean)/std
print("z1 = " ,z1)
print("z2 = " ,z2)
print("0.3621*2=",0.3621*2)
print("For the year 2012, 72,4% of bike rentals are between 107369 and 147369 per month")
Regarding thoses results. I decided to keep the following features :
I did not keep the feature "registered" because it has a correlation of 1 with the "ctn" I did not keet the feature "atemp" because it has a correlation of 1 with the "temp"
ml = file[['season','hr','mnth','yr','weekday','workingday','weathersit','windspeed','temp','cnt']]
X = ml.drop(['cnt'], axis=1)
y = ml['cnt']
from sklearn.model_selection import train_test_split
from sklearn.datasets import load_boston
from sklearn.model_selection import GridSearchCV
from sklearn.metrics import mean_absolute_error
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.ensemble import RandomForestRegressor
from xgboost import XGBRegressor
from vecstack import stacking
from sklearn import ensemble
from sklearn.metrics import r2_score
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)
J'ai choisi d'utiliser XGB. Par conséquent, avec un système d'arbre, je n'ai pas besoin de normaliser les données.
parameters = {'nthread':[4], #when use hyperthread, xgboost may become slower
'objective':['reg:squarederror'],
'learning_rate': [.01,.03, 0.05, .07, .1], #so called `eta` value
'max_depth': [5, 6, 7, 9],
'min_child_weight': [4],
'silent': [1],
'subsample': [0.7],
'colsample_bytree': [0.7],
'n_estimators': [100,500,800,1000]}
xgb1 = XGBRegressor()
xgb_grid = GridSearchCV(xgb1,
parameters,
cv = 2,
n_jobs = -1,
verbose=True)
xgb_grid.fit(X_train,y_train)
print(xgb_grid.best_score_)
print(xgb_grid.best_params_)
model = XGBRegressor(colsample_bytree= 0.7, learning_rate= 0.01, max_depth= 9, min_child_weight= 4, n_estimators=1000, nthread= 4, objective= 'reg:squarederror', silent=1, subsample= 0.7)
model.fit(X_train, y_train)
print(model.score(X_test,y_test))
print(model.score(X_train,y_train))
from sklearn.metrics import mean_squared_error
ypred = model.predict(X_test)
mse = mean_squared_error(y_test,ypred)
print("MSE: %.2f" % mse)
print("RMSE: %.2f" % np.sqrt(mse))
from xgboost import plot_importance
plt.rcParams["figure.figsize"] = (14, 14)
plot_importance(model)
from numpy import sort
from sklearn.feature_selection import SelectFromModel
from sklearn.metrics import accuracy_score
y_pred = model.predict(X_test)
accuracy = r2_score(y_test, y_pred)
print("Accuracy: %.2f%%" % (accuracy * 100.0))
# Fit model using each importance as a threshold
thresholds = sort(model.feature_importances_)
for thresh in thresholds:
# select features using threshold
selection = SelectFromModel(model, threshold=thresh, prefit=True)
select_X_train = selection.transform(X_train)
# train model
selection_model = XGBRegressor(colsample_bytree= 0.7, learning_rate= 0.01, max_depth= 9, min_child_weight= 4, n_estimators=1000, nthread= 4, objective= 'reg:squarederror', silent=1, subsample= 0.7)
selection_model.fit(select_X_train, y_train)
# eval model
select_X_test = selection.transform(X_test)
y_pred = selection_model.predict(select_X_test)
predictions = [round(value) for value in y_pred]
accuracy = r2_score(y_test, predictions)
print("Thresh=%.3f, n=%d, Accuracy: %.2f%%" % (thresh, select_X_train.shape[1], accuracy*100.0))
X = X.drop(['yr','workingday','season'],axis=1)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)
xgb1 = XGBRegressor()
xgb_grid = GridSearchCV(xgb1,
parameters,
cv = 2,
n_jobs = -1,
verbose=True)
xgb_grid.fit(X_train,y_train)
print(xgb_grid.best_score_)
print(xgb_grid.best_params_)
model = XGBRegressor(colsample_bytree= 0.7, learning_rate= 0.03, max_depth= 6, min_child_weight= 4, n_estimators=800, nthread= 4, objective= 'reg:squarederror', silent=1, subsample= 0.7)
model.fit(X_train, y_train)
print(model.score(X_test,y_test))
print(model.score(X_train,y_train))
from sklearn.metrics import mean_squared_error
ypred = model.predict(X_test)
mse = mean_squared_error(y_test,ypred)
print("MSE: %.2f" % mse)
print("RMSE: %.2f" % np.sqrt(mse))
hr = []
mnth= []
weekday= []
weathersit= []
windspeed= []
temp= []
for i in range(0,24):
hr.append(i)
mnth.append(6)
weekday.append(6)
weathersit.append(1)
windspeed.append(0.2)
temp.append(0.4)
var = pd.DataFrame()
var['hr'] = hr
var['mnth'] = mnth
var['weekday'] = weekday
var['weathersit'] = weathersit
var['windspeed'] = windspeed
var['temp'] = temp
var
prediction = model.predict(var)
print(prediction)
fig = go.Figure()
fig.update_layout(
title="Prediction de la location des vélos par heure",
xaxis_title="Heure",
yaxis_title="Nbr de locations",
)
fig.add_trace(go.Scatter(
x=var['hr'],
y=prediction,
name="Prediction"
))
fig.show()
J'avais dans l'idée d'utiliser le stacking afin de créer mon modèle de ML. Cependant, je me suis rendu compte qu'un seul modèle de boosting en choisissant correctement les features au préalable et en éléminant les outliers suffisait. Le boosting permet de diminuer le biais du modèle. Quant au bagging, qui diminue la variance, je m'étais déjà occupé de supprimer les outliers dans le traitement des données.
Pour trouver les meilleurs hyper-parameters, j'ai utilisé GridSearchCV afin d'appliquer de la cross validation. De plus, je me suis penché sur les corrélations entre la feature à prédire "cnt" et les autres afin d'empecher l'overfitting avec "casual" et "registered". Pour finir, j'ai utilisé la feature importance afin d'optimiser mon modèle au mieux. Aussi bien d'un point de vue vitesse d'execution que précision.
Le but premier était de trouver les patterns du dataset, les corrélations, et supprimer les données abérrantes. Comprendre les données. Tout cela pour entreprendre la création d'un modèle dans les meilleures conditions.
Il pourrait être intéressant de comprendre quel genre de personne s'abonne au service et qui préfère rester un simple utilisateur. Afin de faire de la segmentation sur la clientèle et cibler les personnes encore rétisantes à prendre un abonnement. Pour cela, on aurait besoin de d'autres features comme l'âge, la location, les trajets ... des personnes.